## This program generates Figure 4 from "Financial Constraints, Sectoral Heterogeneity,
##  and the Cyclicality of Investment" by Cooper Howes

##### Step 0: Preliminaries #####

## Clear variables and data
rm(list=ls())

## Load packages
library(sandwich)
library(DataCombine)
library(MASS)
library(pracma)
library(mFilter)
library(tidyverse)
library(dynlm)
library(seasonal)
library(glue)

## Set maximum length of response horizon in quarters (default:16)
numres <- 16 

## Set number of standard deviations to show in errors (default: 1.65, i.e. 90% confidence intervals)
se_scale <- 1.65

## Set scale of shock in percentage points (default: 1)
scale <- 1 

## Set number of lags to use in regressions
lags <- 8

## Specify shocks to use in regression (default: "rr", or Romer & Romer (2004)-style narrative shocks)
shocktype <- "rr"

## Regression includes linear time trend, seasonal term, and lags of DV, shock, UR, and GDP growth
f <- formula( glue(paste("L(y,-i) ~ season(y) + {shocktype}_shock",
"+ L({shocktype}_shock,1:lags)  ",
" +L(y,1:lags) +L(ur,1:lags) +L(gdpgrowth,1:lags)  " )))


## Set beginning and ending dates for regressions in (yyyy,q) format
begin <- c(1970,1) ## Default: c(1970,1)
fin   <- c(2008,4) ## Default: c(2008,4)


## Load data from 1966-2012
data <- ts(read.csv("Aggregate_data.csv"),start=c(1966,1),end=c(2021,4),frequency=4)


##### Step 2: Calculate user cost and its components #####


## Create measures to be used in regressions
gdpgrowth <- log(data[,"GDP"]) - stats::lag(log(data[,"GDP"]),-1)
inflation <- log(data[,"gdppi"]) - stats::lag(log(data[,"gdppi"]),-1)
invest_inflation <- log(data[,"Nonres_fixed_invest_price"]) - stats::lag(log(data[,"Nonres_fixed_invest_price"]),-1)
relprice <- data[,"Nonres_fixed_invest_price"]/data[,"ppi"]

logprice <- log(relprice)

avginflation = 4*mean(window(inflation,start=c(1970,1),end=c(2008,4)))
avginvestinflation = mean(window(invest_inflation,start=c(1970,1),end=c(2008,4)))
avgdeprec = 4*mean(window(data[,"bea_deprec_all"],start=c(1970,1),end=c(2008,4)))

ir1 <- 1+data[,"ffr"]/100
ir2 <- 1+data[,"baa_yield"]/100
ir3 <- 1+data[,"adjusted_compustat_total_ir"]

## Adjust for the highest statutory tax rate
tax <- data[,"corptax"]/100


## Weighted financing costs taken from Chirinko et al
## Note that "adjusted_compustat_total_ir" is from Compustat and excluded from the replication package
fcost <- 0.67*(data[,"divyield"]/100 + 0.024) + 0.33*(data[,"adjusted_compustat_total_ir"]*(1-tax) - avginflation)


uc <- relprice*(data[,"bea_deprec_all"] + fcost - avginvestinflation)


## Specify which variables to use in regressions
outcomenames <- c("ir1","fcost","relprice","uc")
longnames <- c("Fed funds rate","Financing cost","Relative investment price","User cost")

## Create empty arrays to hold coefficients and standard errors
gamma    <- matrix(0,nrow=(numres+1),ncol=length(outcomenames))
gamma_se <- matrix(0,nrow=(numres+1),ncol=length(outcomenames))



##### Step 3: Estimate and plot IRFs #####


## Create new window
windows(width=25,height=15)
par(oma=c(5,2.5,0,0))
m <- matrix(1:length(outcomenames),nrow = 2,byrow=T)
layout(mat = m)


## Outer loop over variables to be estimated
for(j in outcomenames){

  ## Interest rates are in levels, everything else in logs
  if(grepl("ir",j)){
    
    eval(parse(text=paste0("y <-",j)))
    
  } else {
    
    eval(parse(text=paste0("y <- log(",j,")")))
    
  }
  

  ## Inner loop over response horizons
  for(i in 0:numres){
    
    ## Need to calculate date of last shock used so all data end in 2008. For
    ##  example, if the data end in 2008, then the latest shock used when estimating
    ##  the response 13 quarters after the shock will need to be in the third
    ##  quarter of 2005, since (2005Q3 + 13 quarters = 2008Q4).
    lastshock_y <- fin[1]-floor(i/4)
    lastshock_q <- 4-i%%4
    
    ## Estimate responses and save coefficients and SEs
    eq <- dynlm(formula=f,start=begin,end=fin,data=data)
    newey <- NeweyWest(eq,prewhite=FALSE)
    
    gamma[(1+i),which(outcomenames==j)] <- coef(summary(eq))[paste0(shocktype,"_shock"),"Estimate"]
    gamma_se[(1+i),which(outcomenames==j)] <- sqrt(newey[paste0(shocktype,"_shock"),paste0(shocktype,"_shock")])
    
  }
  
  
  ## Create IRFs and SEs from array of coefficient estimates
  irf <- gamma[,which(outcomenames==j)]
  irf_se <- gamma_se[,which(outcomenames==j)]
  
  irf1 <- irf
  irf2 <- irf + se_scale*irf_se
  irf3 <- irf - se_scale*irf_se
  
  ## Automatically generate axis ranges
  ymin<-min(irf1,irf2,irf3)
  ymax<-max(irf1,irf2,irf3)
  y_range <- c(ymin,ymax)
  
  ## Generate x-axis sequence for plots
  qtrs <- c(0:numres) 
  
  ## Plot coefficient estimates and SE bands
  par(mar = c(2,2.5,1.5,2.5))
  plot(x=qtrs,y=irf1,type="l",col="blue",lwd=3,ylim=y_range,ann=FALSE,cex.main=1.25,
       axes=F,frame=T)
  axis(1,cex.axis=1.5)
  
  ## Separately label axes for interest rate (pp) and other variables (%)
  if(j=="ir1") {
    axis(2,cex.axis=1.5,at=pretty(y_range),lab=paste0(pretty(y_range*100), "pp"),las=TRUE)} else {
      axis(2,cex.axis=1.5,at=pretty(y_range),lab=paste0(pretty(y_range*100), "%"),las=TRUE)
    }
  
  lines(x=qtrs,y=irf2,type="l",lty=2,col="red",lwd=2)
  lines(x=qtrs,y=irf3,type="l",lty=2,col="red",lwd=2)
  abline(a=0,b=0,col="black",lwd=1)
  if(j %in% outcomenames[3:4]){mtext(text="Quarters after shock",side=1,line=2.5,cex=1.25)}
  title(main=longnames[which(outcomenames==j)], col.main="blue", font.main=2,cex.main=1.5)
  
  
}

## Add legend
par(fig = c(0, 1, 0, 1), oma = c(0, 0, 0, 0), mar = c(0, 0, 0, 0), new = TRUE)
plot(0, 0, type = "n", bty = "n", xaxt = "n", yaxt = "n")
legend(x="bottom",inset=0,legend=c("Coefficient estimate","90% confidence interval"),
       lwd=c(2,2),lty=c(1,2),col=c("blue","red"),cex=1.5,horiz=TRUE)

